########################################################
# Replication file for Ciudad Autonoma de Buenos Aires #
########################################################


rm(list=ls(all=TRUE))
#setwd("")
library(eiPack)
library(coda)
library(foreign)

data=read.csv(file.path(getwd(),"data","CABA_data.csv"))


# Define row and column marginals of the transition matrix. 
p3007=data$paso_3007
p502o=data$paso_502-data$paso_3007
p501=data$paso_501
p503=data$paso_503
pother=data$paso_8+data$paso_68+data$paso_187+data$paso_234+data$paso_324+data$paso_325+data$paso_505+data$paso_506+data$paso_262+data$paso_504
pno=data$votantes-p3007-p502o-p501-p503-pother

g501=data$gral_501
g502=data$gral_502
g503=data$gral_503
gother=data$gral_187+data$gral_505+data$gral_506
gno=data$votantes-g501-g502-g503-gother


data1=data.frame(p3007,p502o,p501,p503,pother,pno,g501,g502,g503,gother,gno)



# Tunning the EI algorithm
tune.nocov=tuneMD(cbind(g501,g502,g503,gother,gno)~cbind(p3007,p502o,p501,p503,pother,pno),
                  data=data1,ntunes=10,totaldraws=200000)


# MCMC
out.nocov=ei.MD.bayes(cbind(g501,g502,g503,gother,gno)~cbind(p3007,p502o,p501,p503,pother,pno),
                      covariate=NULL,data=data1,tune.list=tune.nocov,verbose=0,ret.beta="r",thin=80,burnin=100000,ret.mcmc=TRUE) 


#save(out.nocov, file = "CABAMCMC.Rdata") # Saves the MCMC runs

# Obtain transition matrix for each mesa: 
source(file.path(getwd(),"Codes","getbetas.R"))
beta11=getbetas(out.nocov,"p3007","g502")
beta12=getbetas(out.nocov,"p3007","g501")
beta13=getbetas(out.nocov,"p3007","g503")
beta14=getbetas(out.nocov,"p3007","gother")
beta15=getbetas(out.nocov,"p3007","gno")
beta21=getbetas(out.nocov,"p502o","g502")
beta22=getbetas(out.nocov,"p502o","g501")
beta23=getbetas(out.nocov,"p502o","g503")
beta24=getbetas(out.nocov,"p502o","gother")
beta25=getbetas(out.nocov,"p502o","gno")
beta31=getbetas(out.nocov,"p501","g502")
beta32=getbetas(out.nocov,"p501","g501")
beta33=getbetas(out.nocov,"p501","g503")
beta34=getbetas(out.nocov,"p501","gother")
beta35=getbetas(out.nocov,"p501","gno")
beta41=getbetas(out.nocov,"p503","g502")
beta42=getbetas(out.nocov,"p503","g501")
beta43=getbetas(out.nocov,"p503","g503")
beta44=getbetas(out.nocov,"p503","gother")
beta45=getbetas(out.nocov,"p503","gno")
beta51=getbetas(out.nocov,"pother","g502")
beta52=getbetas(out.nocov,"pother","g501")
beta53=getbetas(out.nocov,"pother","g503")
beta54=getbetas(out.nocov,"pother","gother")
beta55=getbetas(out.nocov,"pother","gno")
beta61=getbetas(out.nocov,"pno","g502")
beta62=getbetas(out.nocov,"pno","g501")
beta63=getbetas(out.nocov,"pno","g503")
beta64=getbetas(out.nocov,"pno","gother")
beta65=getbetas(out.nocov,"pno","gno")


betaout=data.frame(beta11,beta12,beta13,beta14,beta15,beta21,beta22,beta23,
                   beta24,beta25,beta31,beta32,beta33,beta34,beta35,beta41,
                   beta42,beta43,beta44,beta45,beta51,beta52,beta53,beta54,
                   beta55,beta61,beta62,beta63,beta64,beta65,data$votantes,
                   data$comuna,data$mesa)
write.table(betaout,file=file.path(getwd(),"Output","CABA_Mesa.txt"),col.names=NA)


# Obtain aggregate transition matrix. 
source(file.path(getwd(),"Codes","getbetasagr.R"))
beta11=getbetasagr(out.nocov,"p3007","g502")
beta12=getbetasagr(out.nocov,"p3007","g501")
beta13=getbetasagr(out.nocov,"p3007","g503")
beta14=getbetasagr(out.nocov,"p3007","gother")
beta15=getbetasagr(out.nocov,"p3007","gno")
beta21=getbetasagr(out.nocov,"p502o","g502")
beta22=getbetasagr(out.nocov,"p502o","g501")
beta23=getbetasagr(out.nocov,"p502o","g503")
beta24=getbetasagr(out.nocov,"p502o","gother")
beta25=getbetasagr(out.nocov,"p502o","gno")
beta31=getbetasagr(out.nocov,"p501","g502")
beta32=getbetasagr(out.nocov,"p501","g501")
beta33=getbetasagr(out.nocov,"p501","g503")
beta34=getbetasagr(out.nocov,"p501","gother")
beta35=getbetasagr(out.nocov,"p501","gno")
beta41=getbetasagr(out.nocov,"p503","g502")
beta42=getbetasagr(out.nocov,"p503","g501")
beta43=getbetasagr(out.nocov,"p503","g503")
beta44=getbetasagr(out.nocov,"p503","gother")
beta45=getbetasagr(out.nocov,"p503","gno")
beta51=getbetasagr(out.nocov,"pother","g502")
beta52=getbetasagr(out.nocov,"pother","g501")
beta53=getbetasagr(out.nocov,"pother","g503")
beta54=getbetasagr(out.nocov,"pother","gother")
beta55=getbetasagr(out.nocov,"pother","gno")
beta61=getbetasagr(out.nocov,"pno","g502")
beta62=getbetasagr(out.nocov,"pno","g501")
beta63=getbetasagr(out.nocov,"pno","g503")
beta64=getbetasagr(out.nocov,"pno","gother")
beta65=getbetasagr(out.nocov,"pno","gno")
betaoutagr=data.frame(beta11,beta12,beta13,beta14,beta15,beta21,beta22,beta23,
                   beta24,beta25,beta31,beta32,beta33,beta34,beta35,beta41,
                   beta42,beta43,beta44,beta45,beta51,beta52,beta53,beta54,
                   beta55,beta61,beta62,beta63,beta64,beta65,data$votantes,
                   data$comuna,data$mesa)

betaoutagr=as.vector(betaoutagr[1,1:120])
betaoutagr=matrix(betaoutagr,4)
betaoutagr_coefs=t(matrix(betaoutagr[1,],5))
betaoutagr_sd=t(matrix(betaoutagr[2,],5))
write.table(betaoutagr_coefs,file=file.path(getwd(),"Output","CABAAggregate_Coefs.txt"),col.names=NA)
write.table(betaoutagr_sd,file=file.path(getwd(),"Output","CABAAggregate_SD.txt"),col.names=NA)



# Obtain transition matrices for each Comuna. 
source(file.path(getwd(),"Codes","getbetascom.R"))
beta11=getbetascom(out.nocov,"p3007","g502")
beta12=getbetascom(out.nocov,"p3007","g501")
beta13=getbetascom(out.nocov,"p3007","g503")
beta14=getbetascom(out.nocov,"p3007","gother")
beta15=getbetascom(out.nocov,"p3007","gno")
beta21=getbetascom(out.nocov,"p502o","g502")
beta22=getbetascom(out.nocov,"p502o","g501")
beta23=getbetascom(out.nocov,"p502o","g503")
beta24=getbetascom(out.nocov,"p502o","gother")
beta25=getbetascom(out.nocov,"p502o","gno")
beta31=getbetascom(out.nocov,"p501","g502")
beta32=getbetascom(out.nocov,"p501","g501")
beta33=getbetascom(out.nocov,"p501","g503")
beta34=getbetascom(out.nocov,"p501","gother")
beta35=getbetascom(out.nocov,"p501","gno")
beta41=getbetascom(out.nocov,"p503","g502")
beta42=getbetascom(out.nocov,"p503","g501")
beta43=getbetascom(out.nocov,"p503","g503")
beta44=getbetascom(out.nocov,"p503","gother")
beta45=getbetascom(out.nocov,"p503","gno")
beta51=getbetascom(out.nocov,"pother","g502")
beta52=getbetascom(out.nocov,"pother","g501")
beta53=getbetascom(out.nocov,"pother","g503")
beta54=getbetascom(out.nocov,"pother","gother")
beta55=getbetascom(out.nocov,"pother","gno")
beta61=getbetascom(out.nocov,"pno","g502")
beta62=getbetascom(out.nocov,"pno","g501")
beta63=getbetascom(out.nocov,"pno","g503")
beta64=getbetascom(out.nocov,"pno","gother")
beta65=getbetascom(out.nocov,"pno","gno")

betaoutcom=rbind(beta11,beta12,beta13,beta14,beta15,beta21,beta22,beta23,
                   beta24,beta25,beta31,beta32,beta33,beta34,beta35,beta41,
                   beta42,beta43,beta44,beta45,beta51,beta52,beta53,beta54,
                   beta55,beta61,beta62,beta63,beta64,beta65)
betaoutcom=array(betaoutcom,dim=c(5,6,15))
write.table(betaoutcom,file=file.path(getwd(),"Output","CABAComuna_coef.txt"),col.names=NA)




# save(out.nocov, file = "MainResultMCMC.Rdata") # Save MCMC runs

# Geweke Convergence Statistics
geweke=geweke.diag(out.nocov$draws$Beta)
gewekestats=geweke[[1]]
gewekedf=data.frame(gewekestats)
#write.table(gewekedf,file="gewekelist.txt",col.names=NA)


jpeg(file.path(getwd(),"Output","CABAHistogramGeweke.jpeg"))
hist(gewekestats)
dev.off()

print("Max Geweke")
print(max(gewekestats))
print("Min Geweke")
print(min(gewekestats))
print("Percentage Above 2.35")
print(length(gewekestats[abs(gewekestats)>2.35])/length(gewekestats))
print("Percentage Above 1.96")
print(length(gewekestats[abs(gewekestats)>1.96])/length(gewekestats))
print("Percentage Above 1.5")
print(length(gewekestats[abs(gewekestats)>1.5])/length(gewekestats))
print("Percentage Above 1.5")
print(length(gewekestats[abs(gewekestats)>1])/length(gewekestats))
print("Percentage Above 0.5")
print(length(gewekestats[abs(gewekestats)>0.5])/length(gewekestats))
